{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1> Polymorphic Epitope Prediction </h1>\n",
    "\n",
    "This tutorial illustrates how to use Fred2 to integrate genomic information of a patient for epitope prediction\n",
    "\n",
    "This tutorial entails:\n",
    "- Variant construction\n",
    "- Polymorphic Transcript/Protein/Peptide construction\n",
    "- Polymorphic epitope prediction\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "%reload_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "from Fred2.Core import Allele, Peptide, Protein,generate_peptides_from_proteins\n",
    "from Fred2.IO import read_lines, read_fasta\n",
    "from Fred2.EpitopePrediction import EpitopePredictorFactory\n",
    "from Fred2.Core import generate_transcripts_from_variants, generate_proteins_from_transcripts \n",
    "from Fred2.Core import generate_peptides_from_variants\n",
    "from Fred2.IO import read_annovar_exonic\n",
    "from Fred2.IO import MartsAdapter\n",
    "from Fred2.IO import EIdentifierTypes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2> Chapter 1: Generating polymorphic epitopes </h2>\n",
    "<br/>\n",
    "We first have to construct variants to work with. We either can do this manually or by using one of the IO functions of FRED2. Currently, FRED2 supports annotated exonic ANNOVAR files. Once the variant objects are created, we can use them to construct polymorphic transcripts with `Fred2.Core.generate_transcript_from_variants`. For that we also have to specify from which database the additional transcript information (like sequence etc.) should be extracted and what kind of identification system (e.g. RefSeq, ENSEMBL) was used to annotate the variants. Here we use the `Fred2.IO.MatrsAdapter` to connect to the remote BioMart DB and use `RefSeq` as indetification system via specifying this with `Fred2.IO.EIdentifierTypes.REFSEQ` We can also specify which of the community BioMart DB should be used instead of the central BioMart server with the named argument `biomart=`.<br/>\n",
    "\n",
    "`Fred2.Core.generate_transcript_from_variants` will generate all combinatorial possibilities of heterozygous and homozygous variants and incorporate these into the reference transcript sequence. This also means that the function becomes quickly unpractical once a large amount of heterozygous variants should be processed. $n$ heterozygous variants will generate $2^n$ transcript objects. This procedure is done, since no phasing information of the heterozygous variants are routinely available.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "vars = read_annovar_exonic(\"./data/test_annovar.out\")\n",
    "mart =  MartsAdapter(biomart=\"http://grch37.ensembl.org/biomart/martservice?query=\")\n",
    "trans = generate_transcripts_from_variants(vars, mart, EIdentifierTypes.REFSEQ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once we generated the polymorphic transcripts, we can use `Fred2.Core.generate_proteins_from_transcripts` to construct protein sequences. The so generated protein sequences will contain the non-synonymous variants that effected its protein sequence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "proteins = generate_proteins_from_transcripts(trans)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By using `Fred2.Core.generate_peptides_from_proteins`, we can now generate polymorphic peptide objects from the previously generated polymorphic proteins. In addition to the proteins from which peptides are be generate, we have to specify a desired peptide length (e.g. 9-mers)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "peptides1 = list(generate_peptides_from_proteins(proteins, 9))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we are only interested in polymorphic peptides, or the high number of heterozygous variants prohibits the construction of all polymorphic transcripts/proteins, we can use `Fred2.Core.generate_peptides_from_variants`. This function restricts the combinatorial exploration to a specific window size, thereby reducing the number of possible combination to $2^m$ with $m << n$. The window size represents the length of the desired peptides (e.g. 9-mers)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "peptides2 = generate_peptides_from_variants(vars, 9, mart, EIdentifierTypes.REFSEQ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "source": [
    "**Note**: All function starting with `generate` are true generator functions. That means, they stall the calculations until the actual objects are needed, but can only be traversed once."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2> Chapter 2: Epitope prediction </h2>\n",
    "<br/>\n",
    "Once we have generated the peptide objects, we can for example predict their binding affinity. For that, we first have to initialize HLA allele objects, and an epitope prediction method. For more information see the <a href=\"https://github.com/FRED-2/Fred2/blob/master/Fred2/tutorials/EpitopePrediction.ipynb\">tutorial on epitope prediction</a>.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>A*02:01</th>\n",
       "      <th>B*15:01</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Seq</th>\n",
       "      <th>Method</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>(A, A, A, Q, E, A, Q, A, D)</th>\n",
       "      <th>svmhc</th>\n",
       "      <td>-0.995096</td>\n",
       "      <td>-1.295513</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(A, A, D, F, P, R, W, K, R)</th>\n",
       "      <th>svmhc</th>\n",
       "      <td>-1.068312</td>\n",
       "      <td>-0.813017</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(A, A, F, L, A, G, L, L, S)</th>\n",
       "      <th>svmhc</th>\n",
       "      <td>-1.558208</td>\n",
       "      <td>-0.920914</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(A, A, F, L, L, Q, H, V, Q)</th>\n",
       "      <th>svmhc</th>\n",
       "      <td>-1.317206</td>\n",
       "      <td>-1.189752</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(A, A, F, M, Y, V, F, Y, V)</th>\n",
       "      <th>svmhc</th>\n",
       "      <td>0.607884</td>\n",
       "      <td>-0.709952</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                     A*02:01   B*15:01\n",
       "Seq                         Method                    \n",
       "(A, A, A, Q, E, A, Q, A, D) svmhc  -0.995096 -1.295513\n",
       "(A, A, D, F, P, R, W, K, R) svmhc  -1.068312 -0.813017\n",
       "(A, A, F, L, A, G, L, L, S) svmhc  -1.558208 -0.920914\n",
       "(A, A, F, L, L, Q, H, V, Q) svmhc  -1.317206 -1.189752\n",
       "(A, A, F, M, Y, V, F, Y, V) svmhc   0.607884 -0.709952"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alleles = read_lines(\"./data/alleles.txt\", in_type=Allele)\n",
    "svmhc = EpitopePredictorFactory(\"svmhc\")\n",
    "pred_df = svmhc.predict(filter(lambda x: \"*\" not in str(x), peptides1), alleles=alleles)\n",
    "pred_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2> Chapter 3: Post-processing </h2>\n",
    "<br/>\n",
    "These polymorphic peptides have functionalities that allow the user to identify the variants that influenced the peptide sequences and locate their positions within the peptide sequence. With `Peptide.get_variants_by_protein()` we obtain a list of variants that influenced the peptide sequence originating from a specific protein. `Peptide.get_variants_by_protein_position()` returns a dictionary of where the key is the relative position of a variant within the peptide sequence that originated from a specific protein and protein position."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "NM_144701:FRED2_1 [Variant(g.67705958G>A)]\n",
      "{1: [Variant(g.67705958G>A)]}  Protein position:  379  Peptide:  FQTGIKRRI\n",
      "\n",
      "NM_022162:FRED2_3 [Variant(g.50756540G>C)]\n",
      "{8: [Variant(g.50756540G>C)]}  Protein position:  899  Peptide:  SLQFLGFWR\n",
      "\n",
      "NM_022162:FRED2_1 [Variant(g.50756540G>C)]\n",
      "{8: [Variant(g.50756540G>C)]}  Protein position:  899  Peptide:  SLQFLGFWR\n"
     ]
    }
   ],
   "source": [
    "poly_peps = filter(lambda x: any(x.get_variants_by_protein(prot.transcript_id) for prot in x.get_all_proteins()) , peptides1)\n",
    "c=0\n",
    "for p in poly_peps:\n",
    "    c+=1\n",
    "    if c>=3:\n",
    "        break\n",
    "    for prot,poss in p.proteinPos.iteritems():\n",
    "        print\n",
    "        print prot, p.get_variants_by_protein(prot)\n",
    "        for pos in poss:\n",
    "            vars = p.get_variants_by_protein_position(prot, pos)\n",
    "            if vars:\n",
    "                print vars,\" Protein position: \",pos,\" Peptide: \",p"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
